# OLS marginal effects function
marginal.effects <- function(model, conf=90, beta.coef, interaction=F, interaction.coef=NA, Z=1){

	beta1 <- coef(summary(model))[beta.coef]
	beta2 <- coef(summary(model))[interaction.coef]

	if(conf==80){conf.coef<-1.280}
	if(conf==85){conf.coef<-1.440}
	if(conf==90){conf.coef<-1.645}
	if(conf==95){conf.coef<-1.96}
	if(conf==99){conf.coef<-2.575}

	if(interaction==T){
		point <- beta1 + beta2 * Z
		point.l <- point - conf.coef*(sqrt( vcov(model)[beta.coef,beta.coef] + (Z^2)*vcov(model)[interaction.coef,interaction.coef] + 2*Z*vcov(model)[beta.coef,interaction.coef]))
		point.u <- point + conf.coef*(sqrt( vcov(model)[beta.coef,beta.coef] + (Z^2)*vcov(model)[interaction.coef,interaction.coef] + 2*Z*vcov(model)[beta.coef,interaction.coef]))
		int1 <- c(point, point.l, point.u)

		point <- beta1
		point.l <- point - conf.coef*(sqrt(vcov(model)[beta.coef,beta.coef]))
		point.u <- point + conf.coef*(sqrt(vcov(model)[beta.coef,beta.coef]))
		int0 <- c(point, point.l, point.u)

		return(rbind(int0, int1))
	}

	else{
		point <- beta1
		point.l <- point - conf.coef*(sqrt(vcov(model)[beta.coef,beta.coef]))
		point.u <- point + conf.coef*(sqrt(vcov(model)[beta.coef,beta.coef]))
		return(c(point, point.l, point.u))
	}
}
